MAGIC: Microlensing Analysis Guided by Intelligent Computation

Together with u₀, parameters t₀ and t_E are the so-called Paczynski parameters (Paczynski 1986). Unlike the parameters that define the binary properties, these Paczynski parameters are constrained by the global properties of the light curve. In many low-q events, especially planetary events, one can accurately constrain the Paczynski parameters by masking out the relatively brief anomalous features (e.g., Gould & Loeb 1992). Of all three Paczynski parameters, the time-dependent ones, namely t₀ and t_E, enter the microlensing phenomenon via the parameter τ (see Equation (5)). For these reasons, we develop a stand-alone machine-learning method to infer the time-dependent parameters, which is dubbed locator.

The basis of locator is the U-Net architecture. U-Net is a fully convolutional neural network that is originally designed to perform semantic segmentation of images (Ronneberger et al. 2015). U-Net takes the original image as the input, and outputs a mask of the same shape. The mask can be interpreted as a pixel-wise classification result, representing, for example, the probability of each pixel belonging to some specific segment. U-Net has been ubiquitously used in biomedical image processing (Siddique et al. 2021) since it was proposed in 2015. It has also seen a growing number of applications in astronomy (e.g., Van Oort et al. 2019; Makinen et al. 2021).

The architecture of our locator closely follows that of (Ronneberger et al. 2015). See their Figure 1 for an illustration. In detail, locator differs in three aspects. First, all convolutional layers are one-dimensional, because locator deals with time series data rather than image data. Second, the pure convolutional blocks are replaced by residual blocks (He et al. 2016) with Parametric Rectified Linear Unit (PReLU, He et al. 2015) activation functions to ease gradient back propagation. Third, we set the channel and kernel sizes to 128 and 7, respectively, to better match the shape of the input data. In practice, the performance of locator is not sensitive to the values of these hyperparameters.

We design locator to identify the light-curve segment that falls within a predefined time window ΔT ≡ [t₀ − kt_E, t₀ + kt_E]. Here k is a parameter that can be tuned to optimize the output performance, and its choice will be discussed in Section 3.3. For each light curve, which is irregularly sampled (see Section 3.2), we first interpolate it with natural cubic splines and extract the values at n_t = 4000 ⁴ equally stepped time stamps ${\{{t}_{j}\}}_{j=1}^{{n}_{t}}$. ⁵ These time stamps are assigned labels depending on whether or not they are inside ΔT,

$$\begin{eqnarray}{r}_{j}=\left\{\begin{array}{ll}1, & {t}_{j}\in {\rm{\Delta }}T\\ 0, & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 6 }$$

At run time, locator takes the time series and attempts to produce a mask ${\{{p}_{j}\}}_{j=1}^{{n}_{t}}$ that best approximates the ground truth labels.

We train locator on simulated light curves that match the key properties of real microlensing light curves (see Section 3.2). The loss function that locator minimizes is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{locator}}={{ \mathcal L }}_{\mathrm{BCE}}+{{ \mathcal L }}_{\mathrm{Dice}}+{{ \mathcal L }}_{\mathrm{reg}}.\end{eqnarray} \tag{ 7 }$$

Here ${{ \mathcal L }}_{\mathrm{BCE}}$ is the binary cross entropy (BCE) loss function

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{BCE}}\equiv {\mathbb{E}}[r\mathrm{ln}p+(1-r)\mathrm{ln}(1-p)],\end{eqnarray} \tag{ 8 }$$

which is the standard loss function for binary classification (Goodfellow et al. 2016). The expectation runs over all time stamps of all light curves in the training set. The second term in Equation (7) is the Dice loss function (Sudre et al. 2017), which is introduced to assuage the class imbalance problem in events with small t_E:

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{Dice}}\equiv {\mathbb{E}}\left[1-\displaystyle \frac{\sum _{j}{r}_{j}{p}_{j}}{\sum _{j}({r}_{j}+{p}_{j})}\right].\end{eqnarray} \tag{ 9 }$$

The third term in Equation (7) is the regularization term to prevent the mask from wiggling up and down more than twice

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{reg}}\equiv {\mathbb{E}}{\left(\sum _{j=1}^{{n}_{t}-1}| {p}_{j+1}-{p}_{j}| -2\right)}^{2}.\end{eqnarray} \tag{ 10 }$$

The expectations in the last two terms are performed on all light curves in the training set.

The output mask from locator contains information about the event peak time, t₀, and the timescale, t_E. Specifically, the centroid of the mask corresponds to t₀ and the total enclosed area is directly related to t_E. In the continuum form, they are given by

$$\begin{eqnarray}&&{t}_{0}=\displaystyle \frac{\int p(t){tdt}}{\int p(t){dt}},\quad {t}_{{\rm{E}}}=\displaystyle \frac{1}{2k}\int p(t){dt}.\end{eqnarray} \tag{ 11 }$$

These integrals are evaluated via the trapezoidal rule in the discrete time series.

To simulate light curves for the training of locator, we randomly sample 2000 time stamps out of a total duration of 270 days. This corresponds to an averaged cadence of ∼3 hr, with the lower and upper 10th percentiles corresponding to 20 minutes and 7.2 hr. For each light curve, we randomly draw t_E from a truncated log-normal distribution with a mean of 10^1.15 days, a standard deviation of 10^0.45 days, and boundaries at 5 and 100 days. This log-normal distribution matches approximately the observed distribution of t_E (e.g., Mróz et al. 2017), and the truncations are introduced to eliminate relatively short-timescale and long-timescale events, both of which may be subject to additional issues (e.g., parallax effect for long events) and are not the focus of the present work. As for t₀, even though it can be any value during (and possibly outside) the sampled time interval, it is reasonable to believe that a very rough by-eye estimate can limit t₀ to within ∼±t_E of the truth. We therefore randomly draw t₀ out of a uniform distribution between − t_E and t_E. To summarize,

$$\begin{eqnarray*}&&\begin{array}{l}{t}_{{\rm{E}}}/\mathrm{day}\sim \mathrm{TruncLogNorm}(5,100,\mu ={10}^{1.15},\sigma ={10}^{0.45}),\\ \ {t}_{0}\sim \mathrm{Uniform}(-{t}_{E},{t}_{E}).\end{array}\end{eqnarray*}$$

Although the goal of locator is to infer t₀ and t_E and not the other parameters, the training data set should be representative of the target—the binary light curves. We therefore sample the remaining parameters in the following way

$$\begin{eqnarray*}\begin{array}{rcl}{u}_{0} & \sim & \mathrm{Uniform}(0,1),\\ q & \sim & \mathrm{LogUniform}({10}^{-3},{10}^{0}),\\ s & \sim & \mathrm{LogUniform}(0.3,3),\\ \alpha & \sim & \mathrm{Uniform}(0,360),\\ {f}_{{\rm{S}}} & \sim & \mathrm{LogUniform}({10}^{-1},{10}^{0}),\end{array}\end{eqnarray*}$$

and ρ is fixed at 10⁻³, given that the finite-source effect is not expected to be detected for the adopted (effective) cadence. Note that the mass ratio q is limited to > 10⁻³ simply because we want to focus on relatively massive binaries. With increased sampling cadence locator can in principle be extended down to lower q values. We defer a detailed exploration to future works.

With the sampled binary parameters and time stamps, we calculate the microlensing light curve using the VBBL algorithm (Bozza 2010; Bozza et al. 2018). To introduce noises to the simulated light curves, we adopt a simplified noise model. Gaussian noises with S/Ns of 33 in flux are assumed for all data points, and the flux values are converted to magnitudes for a baseline object with m_base = 18.

Simulated light curves without significant binary features, defined as ${\chi }_{\mathrm{pspl}}^{2}$/d.o.f. <1.25, ⁶ are rejected. In the end, two data sets each containing 10⁵ light curves are generated, one for training and the other for validation and test.

For our default run, we set k = 1/3 so that ΔT with the longest timescale (100 days) can still be fully covered in the simulation interval [−135, 135]. We train locator with the ADAM optimizer (Kingma & Ba 2014). The learning rate is set to 4 × 10⁻³ initially and drops by ten percent after each epoch. Training samples are divided into minibatches of size 128. After 13 epochs (∼2 hr) of training on one NVIDIA Tesla V100 GPU, the average loss drops to 0.16 on the validation set of size 1024. The optimized model is then used to predict the mask and estimate t₀ and t_E according to Equation (11).

We evaluate locator on 16,384 light curves from the test set. The output masks are illustrated in Figures 2 and 3, and the comparisons between the predicted and ground truth values of lgt_E and t₀/t_E are shown in Figure 4. The accuracy of the prediction is measured by the rms error (RMSE) between the ground truth and the prediction. They are indicated in the upper left corner respectively. We note that locator has very different performances on close (s < 1) binaries and wide (s > 1) binaries. This is especially true for the t₀/t_E parameter. We take the results of close binaries as a more reliable evaluation of the overall performance and defer to the end of this section a more detailed explanation of this phenomenon.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the default run, locator achieves an accuracy in lgt_E with RMSE = 0.049, which translates to a fractional uncertainty of 11.3% on t_E. The accuracy of the ratio t₀/t_E is comparable at RMSE (s < 1) = 0.062, which corresponds to an uncertainty of 1.9 days on t₀ for a typical event with t_E = 30 days. An example event is shown in Figure 5, in which the light curve substantially deviates from the Paczynski (1986) model due to the caustic-crossing feature and yet locator successfully recovers t₀ and t_E at high accuracy.

Zoom In Zoom Out Reset image size

As shown in the left panel of Figure 4, the scattering in lgt_E increases for shorter events. This is probably because, as the timescale decreases, fewer data points are available to locator in the parameter inference. For example, with k = 1/3 there are only ∼25 time stamps that fall within ΔT for the shortest (t_E = 5 days) event. To improve the performance of locator on short events, we select light curves with t_E < 40 days and train another model with k = 2. After eight epochs of training with the same learning rate schedule, the loss drops to 0.083 on the validation set. We evaluate the new model on 16,384 light curves with t_E < 40 days. ⁷ The results are shown in Figure 6. The prediction accuracy on the timescale increases by a factor of ∼2, reaching a fractional uncertainty on t_E of 6.2%. The performance on the ratio t₀/t_E remains largely unchanged.

Zoom In Zoom Out Reset image size

The above investigation has two implications. First, it outlines a general strategy to obtain a more accurate t_E. For example, one may divide the overall range of t_E into multiple smaller intervals and train separate locator models. The preliminary estimate of t_E from the default run of locator can then be "polished" by the locator model that is optimized for the narrower t_E interval including the preliminary estimate. Additionally, the comparison in locator performance between two k values also sheds light on how locator works. In the single-lens case, the most constraining power on t₀ comes from the central region with a width of ∼u₀ t_E, whereas the event timescale t_E is best constrained by data points extending into the light-curve wings. This seems to apply to locator as well. As shown in Figure 7, the timescale parameter can be better constrained when a wider mask window (i.e., a larger k) is used, whereas the event peak time t₀ is largely unaffected by the choice of k. More works are needed to better understand the mechanism of locator (and U-Net in general).

Zoom In Zoom Out Reset image size

We have fixed the magnitude of the baseline object at 18 and used S/N = 33 for all data points in all simulated events. In reality, the baseline magnitude varies from event to event, and the photometric uncertainty varies from data point to data point. However, as discussed previously, the choice of the baseline magnitude does not really matter to the performance of locator, as long as the photometric observations cover a long enough time baseline. We therefore only focus on the impact of varying S/N on the performance of locator. Rather than simulating the realistic noise conditions, ⁸ we study the impact of different training strategies. Several data sets with different S/N values are generated following the method in Section 3.2. We then train locator models on light curves with a fixed S/N value and then evaluate the performance on light curves from the same data set. The result of this training strategy is shown as the blue points in Figure 8. Here we only show the RMSE values of t₀ as the trend is very similar for t_E. Next, we train one locator model on events with the lowest quality (S/N = 10) and apply it to all data sets, the result of which is shown in orange in Figure 8. In the third test, we train locator on events with the highest quality (S/N = 1000) and apply it to all data sets. This result is shown in green in the same figure. As expected, locator performs the best if it is trained and evaluated on light curves with the same (or very similar) data quality. However, the performance of locator is only slightly worse if it is trained on the noisiest data set. This is probably because locator, similar to other convolutional neural networks (Goodfellow et al. 2016), is effectively composed of a series of filters learned from the data. When faced with noisy data, it learns to denoise through training, so the relatively large noises in the training data make locator more robust. For the same reason, the accuracy of the predicted t₀ is substantially worse on noisy light curves if locator only saw during the training light curves of the highest quality. For applications to real events, it seems that training one locator on relatively noisy data may be a good strategy that balances efficiency and accuracy.

Zoom In Zoom Out Reset image size

Finally, we discuss the issue related to wide (s > 1) binaries. As pointed out earlier, the match between the ground truth labels and the predictions of locator for wide binaries is not as good as for close binaries. This larger mismatch, however, does not mean that the inference power of locator is reduced for wide binaries. After visually inspecting a large number of wide binary light curves, we find that the mismatch is related to the use of the reference system. When generating the light curves, we followed the common convention and adopted the center of mass of the binary as the origin of the coordinate system. However, for wide binaries, the external shear due to a component at wide separation is better modeled in the reference frame defined on the center of magnification. For s > 1, the offset between the center of mass and the closer center of magnification is (e.g., Di Stefano & Mao 1996; Chung et al. 2005; Penny 2014)

$$\begin{eqnarray}{\rm{\Delta }}x=\displaystyle \frac{1}{1+q}\left(s-\displaystyle \frac{1}{s}\right)\cdot \left\{\begin{array}{ll}q, & \mathrm{if}\,\mathrm{primary}\\ 1, & \mathrm{if}\,\mathrm{secondary}\end{array}\right..\end{eqnarray} \tag{ 12 }$$

For very wide binaries, the event timescale given by the mass of the closest lens also takes over the timescale given by the total mass of the binary. Both can lead to large offsets between the ground truth window and the predicted window. The most extreme example that has the largest deviation in t₀ and t_E with q = 0.99 and s = 2.8 is shown in Figure 9. Once the reference frame defined on the center of magnification is adopted, the match between the ground truth and the prediction of locator is substantially improved (see also Zhang et al. 2021).

Zoom In Zoom Out Reset image size

We develop a separate method to infer the time-independent microlensing parameters, which we reparameterize as

$$\begin{eqnarray}&&\omega \equiv ({u}_{0},\mathrm{lg}q,\mathrm{lg}s,\alpha /180^\circ ,\mathrm{lg}{f}_{{\rm{S}}}).\end{eqnarray} \tag{ 13 }$$

As shown in Figure 1, this estimator consists of two main parts. A neural controlled differential equation (neural CDE) is used to extract features from the microlensing light curve, which is irregularly sampled and may contain large data gaps. The second part, the mixture density network (MDN), infers the probability distribution of the microlensing parameters.

To improve the training speed and reduce the memory requirement, we first apply the log-signature transform (Morrill et al. 2021). This transform outputs a shorter, steadier, and higher dimensional sequence that summarizes the substep information. For the simulated light curves with 500 time stamps (see Section 4.2), we set the log-signature depth n_D = 3 and the shortened length l = 100. The output thus has v = 5 dimensions.

Neural CDE approximates the underlying process of the time series as a differential equation (Kidger et al. 2020; Kidger 2022)

$$\begin{eqnarray}&&{z}_{t}={\xi }_{\theta }({X}_{{t}_{0}})+{\int }_{{t}_{0}}^{t}{f}_{\theta }\left({z}_{t^{\prime} }\right){{dX}}_{t^{\prime} }.\end{eqnarray} \tag{ 14 }$$

Here the latent state z_t is defined on the time interval of the input signal, (t₀, t_n ), with w = 32 dimensions. The continuous signal X_t , defined on the same time interval, is obtained by applying natural cubic splines to the input, discrete signal. ⁹ The initial value and signal derivative of the latent state, ξ_θ and f_θ , are two different neural networks, where the subscript θ denotes the dependence on the learnable machine-learning parameters. By Equation (14), neural CDE gradually extracts features from the signal X_t according to a learned policy f_θ , records its knowledge in the latent state z_t , and outputs a terminal value ${z}_{{t}_{n}}$ that summarizes the useful information of the input signal. We refer interested readers to Kidger et al. (2020) and Kidger (2022) for further details on neural CDE.

The summary features from neural CDE are then fed into the MDN (Bishop & Nasrabadi 2006) to estimate the probability distribution of the microlensing parameters

$$\begin{eqnarray}&&p(\omega | {z}_{{t}_{n}})=\sum _{i=1}^{{n}_{{\rm{G}}}}{\pi }_{i}({z}_{{t}_{n}}){\phi }_{{\bar{\omega }}_{i}({z}_{{t}_{n}}),{{\rm{\Sigma }}}_{i}({z}_{{t}_{n}})}(\omega ).\end{eqnarray} \tag{ 15 }$$

Here ${\phi }_{\bar{\omega },{\rm{\Sigma }}}$ denotes the density of a multivariate Gaussian with mean $\bar{\omega }$ and covariance matrix Σ. The normalized weight π_i , mean ${\bar{\omega }}_{i}$, and covariance matrix Σ_i of each of the n_G multivariate Gaussians are given by the neural network of MDN. Compared to the masked autoregressive flow method of Zhang et al. (2021), MDN provides a straightforward and efficient alternative with explicit density estimation, which is useful for the identification of degenerate solutions and further optimization (see below). We adopt n_G = 12 Gaussians with diagonal covariance matrices. This preserves universality as long as n_G is large enough (Bishop & Nasrabadi 2006).

In the joint framework, estimator, the initial value network ξ_θ , the evolution network f_θ , and MDN are constructed as residual networks (ResNet; He et al. 2016) with three residual blocks and each block consisting of two fully connected layers of width 1024. The evolution network has an additional activation of a hyperbolic tangent to rectify the dynamic. To train the model, we minimize the negative log-likelihood (NLL) loss function

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{estimator}}=\mathop{{\mathbb{E}}}\limits_{\mathrm{light}\ \mathrm{curve}\ }[-\mathrm{ln}p({\omega }_{\mathrm{true}}| {z}_{{t}_{n}})].\end{eqnarray} \tag{ 16 }$$

To further optimize the microlensing parameters, we start from the mean positions of the multivariate Gaussians and minimize the light-curve χ² values via the Nelder & Mead (1965) simplex algorithm. In this step, the source flux parameter f_S is derived analytically for a given set of (q, s, u₀, α) by maximizing the model likelihood (e.g., Poleski & Yee 2019).

For the purpose of learning the time-independent parameters, we generate light curves as functions of the parameter τ. This is equivalent to setting t₀ = 0 and t_E = 1 day. The other parameters are sampled in the same way as in Section 3.2.

Each light curve contains 500 data points randomly drawn within the interval − 2 ≤ τ ≤ 2. For a typical timescale t_E = 20 days, this corresponds to an averaged cadence of ∼4 hr with a range from 24 minutes (lower 10%) to 8.8 hr (upper 10%). To mimic the missing data points due to bad weather and/or instrumental failures, we throw away data points within one randomly chosen light-curve segment out of the 25 equal-duration segments that constitute the whole light curve. This leaves a data gap of 3.2 days for the aforementioned timescale. Light curves and noises are generated in the same way as in Section 3.2. Finally, we only keep light curves with significant binary features, defined as ${\chi }_{\mathrm{pspl}}^{2}$/d.o.f. > 2.

Six data sets are generated, each having 10⁵ light curves. We use five of them for training estimator and the last one for validation and testing. To further understand the performance of estimator on light curves with different cadences, we have also generated data sets with various numbers of data points within the interval −2 ≤ τ ≤ 2. For these additional data sets, we have assumed irregular sampling but no gap.

We train estimator with the ADAM optimizer (Kingma & Ba 2014). The learning rate, initially set to 10⁻⁴, drops by ten percent after each epoch and reaches a minimum of 10⁻⁶. Training samples are labeled by their input parameters ω and divided into minibatches of size 128. The training took ∼1 day for 50 epochs on one NVIDIA Tesla V100 GPU, and the average loss drops to −8.4 on the validation set of size 1024.

We apply the optimized estimator to 16,384 light curves from the test set and obtain the probability distribution of the microlensing parameters given by the summation of 12 multivariate Gaussians. We then apply the Nelder & Mead (1965) simplex algorithm to get the refined microlensing parameters. We show in Figure 10 the cumulative distributions of the model χ² values of two selected sets of parameters, one that has the smallest χ² value and the other that is closest to the input position in the parameter space. The median χ² values are ∼500 for both solutions, comparable to the number of data points in the simulated light curve and thus approximately the degrees of freedom. This suggests that the found solutions can indeed describe the data very well. Even though in ∼20%–30% of cases the model χ² values are substantially large compared to the number of data points, the returned parameters are still good representatives of the input parameters, as will be discussed later. The Appendix shows an example event whose model prediction gives large χ² and yet reasonably accurate binary parameters.

Zoom In Zoom Out Reset image size

We have included a mixture of n_G = 12 Gaussians in the density estimate. Is it flexible enough? To answer this question, we first look at the distribution of the Gaussian weights given by estimator. As shown in Figure 11, the probability distribution is dominated by one Gaussian in ∼50% of the tested events and by the top three Gaussians with the highest weights in nearly all tested events. As for the optimized solutions, 74% and 85% of the minimum χ² and minimum RMSE solutions are given by the top four Gaussians, respectively. Only in ≲1% of the tested events does the Gaussian with the lowest weight lead to either solution. Therefore, we conclude that 12 Gaussians can satisfactorily describe the probability distribution of microlensing model parameters.

Zoom In Zoom Out Reset image size

We compare the estimator predictions before and after the simplex optimization with the input parameters in Figures 12 and 13. For each tested event, we choose as the model prediction the set of parameters that is closest to the input position in the parameter space, in order to minimize the impact of degenerate solutions. The results change only marginally if we pick the Gaussian with the highest weight as the model prediction, as seen in Figure 14. We use the median absolute deviation (MAD) to evaluate the agreement between input and prediction. Compared to RMSE, MAD is more robust against outliers. As shown in Figure 12, the direct output of estimator already contains a set of parameters that is fairly close to the input truth. Among all parameters, binary mass ratio q and separation s are of special interest to the physical interpretation of the event. The direct output of estimator yields a typical deviation of 0.057 (0.0155) for lgq (lgs), corresponding to a fractional uncertainty of 13% (3.6%) on q (s). The accuracy of the prediction is further improved after the optimization, as shown in Figure 13. Specifically, the optimized prediction can achieve median deviations of ∼0.03 and ∼0.006 for lgq and lgs, corresponding to fractional errors of ∼7% and ∼1.4% in q and s, respectively. Overall, the agreement between the input and predicted values is very well, demonstrating the high accuracy of our method on light curves with realistic sampling conditions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have also applied the optimized estimator to light curves in the additional data sets, which have different averaged cadences but no data gap. ¹⁰ The results are shown in the lower right panel of Figure 13 for the optimized prediction. There are a few interesting trends to report. First, the performance of estimator is better, although fairly marginally, on events with 500 points and no data gap. Second, the prediction accuracy decreases with decreasing number of data points as expected, but the degradation is typically small. Specifically, with 200 data points, the prediction accuracy is only worse by ∼10^0.13 − 1 ≈ 35% compared to the default case with 500 data points. The only exception is seen in the case with 100 data points, for which the performance is reduced by a factor of ≳2. This is probably because the averaged cadence becomes too low to constrain binaries with mass ratios down to 10⁻³. ¹¹ Nevertheless, the overall uncertainty in the case of 100 points remains fairly small and acceptable. These results may provide useful guidance to the future application of our method to real light curves with variant cadences and/or the detection of low-mass-ratio planets.

We show in Figure 15 an example event whose anomalous feature is only partially covered by data points. The probability distribution from estimator, represented by the multivariate Gaussians, is already in the neighborhood of the ground truth position. The optimization with Nelder & Mead (1965) algorithm leads to an even better match in microlensing parameters as well as the model of the light curve. Events like this one, which has irregular samplings and large data gaps, pose challenges to machine-learning methods that are based on the standard RNN approach.

Zoom In Zoom Out Reset image size

Our method also identifies multiple solutions in events that suffer model degeneracy. We defer to Section 5.2 for further discussions on this issue.

Combining the methods described in Sections 3 and 4, the joint pipeline of MAGIC works as follows. First, we feed the light curve into locator, which outputs the predicted t₀ and t_E. If t_E is small, we can further refine it using locator with a larger k. Then we use the predicted t₀ and t_E to shift and rescale the light curve. The resulting light curve is processed by estimator, which in return gives the positions and shapes of n_G = 12 Gaussians in the parameter space. If needed, these Gaussians can be further optimized to obtain more accurate predictions (see Section 4.1).

The performance of the combined pipeline is largely determined by the propagation of inaccurate parameters from locator to estimator. We therefore test the joint pipeline on 16,384 light curves from the test set of locator. After applying locator, we transform the light curves using the predicted values of t₀ and t_E. The transformed light curves are then sent to estimator to infer the time-independent microlensing parameters. The set of parameters closest to the position of the input values in the parameter space is adopted as the final prediction.

The comparisons between the input and the predicted values of the time-independent parameters are shown in Figure 16. The typical, fractional deviation in q and s are 25% and 9%, respectively. These deviations remain small, even though they are worse than the estimator-only case by a factor of ∼2. Overall, the joint pipeline achieves fairly high prediction accuracy in the time-independent parameters even with the use of the inaccurate time-dependent parameters. We expect the accuracy to be further improved once the predictions are optimized via the Nelder & Mead (1965) algorithm, as seen in the estimator-only case.

Zoom In Zoom Out Reset image size

A substantial fraction of binary events can be well described by a Paczynski (1986) model plus relatively short binary perturbation. For such events, one may mask out the anomalous feature and fit a single-lens model to obtain precise and accurate peak time t₀ and event timescale t_E. Combining these values with estimator, we can achieve nearly as accurate time-independent parameters as we do in the estimator-only case.

It is worth noting that the light curves generated for locator are similar to but not the same as those used for training estimator. First, locator uses binary light curves with ${\chi }_{\mathrm{pspl}}^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}.\,\gt \,1.25$, whereas estimator uses binaries with ${\chi }_{\mathrm{pspl}}^{2}/{\rm{d}}.{\rm{o}}.{\rm{f}}.\,\gt \,2$. Both correspond to a Δχ² > 500 between single and binary models, under the reasonable assumption that the binary signals are usually local features. ¹² Furthermore, the time window is fixed and the timescale parameters vary in the data set for locator. This means the number of data points within −2 ≤ τ ≤ 2 varies too. Nevertheless, we have shown in Section 4.3 that the performance of estimator is only mildly affected by the various length of input time series.

In addition to its high accuracy, our method demonstrates high efficiency in parameter estimation. Specifically, locator and estimator can process ∼1000 and ∼400 light curves per second on a single GPU, respectively. As shown in Figure 13, the output of estimator is already accurate enough for purposes like statistical studies. If needed, the optimization step can be performed. Compared to the two machine-learning steps, this traditional optimization step is slower, but there are only up to 12 initial positions to be optimized, which altogether may require up to a few thousand light-curve evaluations. For comparison, the traditional approach usually requires optimizations at individual grid points on a fine grid in the (lgq, lgs, α) space (e.g., Albrow et al. 2000; Wang et al. 2022). Our approach therefore represents a speedup by factors of ∼10²–10⁵.

We apply MAGIC to KMT-2019-BLG-0371 (hereafter KB190371 for short), which is a short-timescale (t_E = 6.5 days) binary event with q ≈ 0.08 or 0.12 (Kim et al. 2021). KB190371 is chosen for this demonstration based on several considerations. First, its caustic-crossing feature is prominent and U-shaped, which is typical for binary events. Second, this event suffers from close/wide degeneracy (Griest & Safizadeh 1998; Dominik 1999) ¹³ , thus enabling additional tests of MAGIC. Furthermore, except for the relatively short timescale, the other microlensing parameters are fairly typical in comparison to those used in our simulation.

The available data sets from Kim et al. (2021) are not aligned to the same magnitude system. In order to align the data sets and then apply MAGIC, we mask out the caustic-crossing feature, fit a Paczynski (1986) model to the data, and use the best-fit flux parameters to align all other data sets to the OGLE one. We note that, although here we have made use of the special property of KB190371 to align the data sets, such a data alignment task is the same as deriving the source color and can be performed without reference to specific microlensing models (Yoo et al. 2004; Gould et al. 2010a), as has been practiced in the analysis of many real microlensing events (e.g., Gould et al. 2010b; Zhu et al. 2017).

To apply our joint pipeline, we first preprocess the data and make the light curve similar to those used in the training procedure. To be more specific, we only include data points falling in the time interval of about 100 days centered at the anomaly, and filter out data points with large (>0.1 mag) error bars. KB190371 was observed with a high cadence, so we bin every five consecutive points into one. Finally, the light curve is shifted to the same base magnitude as those used in locator.

Then we feed the preprocessed light curve into locator. Because the timescale is relatively small, we set k = 2. The output mask is rectified into binary values {0, 1} with a threshold of 0.5. This step gives the prediction ${t}_{0}\,(\mathrm{HJD}^{\prime} )\,=8592.391$ and t_E = 6.68 days, both of which are close to the values are given in Kim et al. (2021). We transform the light curve with the predicted t₀ and t_E and only keep the segment within −2 ≤ τ ≤ 2.

The rest of the application follows closely the procedure given in Section 5.1. With the positions of the 12 multivariate Gaussians given by estimator, we perform Nelder & Mead (1965) optimization to further refine the microlensing parameters. Because the light curve shows a prominent finite-source effect, we also include in the optimization the scaled source size ρ.

The results are shown in Figure 17. The top two Gaussians with the highest weights match closely to the close (s < 1) and wide (s > 1) solutions, respectively. The corresponding model light curves from the optimized solutions both match the data reasonably well, with almost indistinguishable χ² values. According to these solutions, the event can be explained by a binary with (lgq, lgs) = (−1.16, −0.057) or (−0.95, 0.180). Given the imperfect alignment as well as the slightly different preprocessing procedures, we consider these values in good match with those given by the discovery paper.

Zoom In Zoom Out Reset image size

Kim et al. (2021) found the source scaled parameter ρ to be (6.46 ± 0.17) × 10⁻³ and (6.99 ± 0.21) × 10⁻³ for the close and wide solutions, respectively. These are quite different from the value (10⁻³) that we have used to train MAGIC. Nevertheless, the optimized parameters yield ρ = 6.6 × 10⁻³ and 7.3 × 10⁻³ for the close and wide solutions, respectively, in good agreement with the values from the discovery paper. In other words, the presence of the finite-source effect in KB190371 does not seem to confuse MAGIC, even though MAGIC has not learned to infer the source scaled parameter ρ. This is probably because the finite-source effect in KB190371 only affects a limited number of data points around the caustic entrance and exit. While the majority of binary microlensing events bear the same feature, some may have their binary signals substantially modified by the finite-source effect (e.g., Hwang et al. 2018). Therefore, it may be necessary to include in MAGIC the ability to infer ρ. See Section 6 for further discussion.

To further demonstrate the capability of MAGIC to handle large data gaps, we exclude observations from the South Africa site of KMTNet. This creates a substantial gap during the exit of the caustic crossing, as shown in the top right panel of Figure 18. The same procedure is applied as before, and the results are shown in Figure 18. The two degenerate solutions are again identified at high accuracy, indicating the robustness of MAGIC in handling large data gaps as well as model degeneracy.

Zoom In Zoom Out Reset image size